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Variational Monte Carlo study of gapless spin liquid in the spin-1/2 XXZ 
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By using the variational Monte Carlo technique, we study the spin-1/2 XXZ antiferromagnetic 
model (with easy-plane anisotropy) on the kagome lattice. A class of Gutzwiller projected fermionic 
states with a spin Jastrow factor is considered to describe either spin liquids (with U{1) or Z 2 
symmetry) or magnetically ordered phases (with q = (0,0) or q = (47r/3,0)). We find that the 
magnetic states are not stable in the thermodynamic limit. Moreover, there is no energy gain 
to break the gauge symmetry from U{1) to Z 2 within the spin-liquid states, as previously found 
in the Heisenberg model. The best variational wave function is therefore the U{1) Dirac state, 
supplemented by the spin Jastrow factor. Furthermore, a vanishing S = 2 spin gap is obtained at 
the variational level, in the whole regime from the XY to the Heisenberg model. 

PACS numbers: 75.10.Jm, 75.10.Kt, 75.40.Mg, 75.50.Ee 


Introduction. Quantum spin liquids with topological 
order and fractional excitations are exotic states of mat¬ 
ter that do not show any local order down to zero tem¬ 
perature.^ Their importance for the field of correlated 
systems is directly related to the connection to uncon¬ 
ventional electron pairing, thus giving a clue to explain 
the mechanism of high-temperature superconductivity.^’^ 
In the last two decades, there have been intensive stud¬ 
ies suggesting that quantum spin liquids might be sta¬ 
bilized at low temperatures in realistic two-dimensional 
frustrated magnetic systems. The spin-1/2 Heisenberg 
antiferromagnetic model on the kagome lattice represents 
one of the most promising examples.^ From the exper¬ 
imental side, the so-called Herbertsmithite shows very 
promising signatures for magnetically disordered phases 
down to extremely low temperatures.^”^ Even more in¬ 
terestingly, many experimental probes suggested the ex¬ 
istence of gapless spin excitations; in particular, neutron 
scattering measurements highlighted the presence of a 
broad continuum of excitations down to small energies. 

For Herbertsmithite, it is widely believed that the gross 
features can be captured by the spin-1/2 Heisenberg an- 
tiferromagnet on the kagome lattice with the nearest- 
neighbor interactions only. This model has been studied 
by several analytical and numerical approaches in recent 
years, with contradicting out comes. In particular, ac¬ 
curate density-matrix renormalization group (DMRG) 
calculations highlighted the possibility that the ground 
state can be a fully gapped Z 2 topological spin liquid. 

A different scenario has been put forward by using vari¬ 
ational Monte Garlo approaches based upon Gutzwiller 
projected fermionic state,which find a gapless spin 
liquid with a competing ground-state energy.^^’^^ 

In order to clarify the nature of the spin-liquid phase, 
several authors considered the effect of different “per¬ 
turbations” to the nearest-neighbor Heisenberg model, 
the most obvious one being a second-neighbor super¬ 
exchange. However, given the lack of the consis¬ 

tent results,it is still not clear if this term helps 


the stabilization of the spin liquid or not. The inclusion of 
an additional third-neighbor couplings^^’^^’^^’^^ or three- 
spin chiral interactions^^ stabilizes a topological spin liq¬ 
uid with spontaneously time-reversal symmetry break- 
ing,35’36 which has been identified as the z/ = 1/2 bosonic 
quantum Hall state. ^^’32 ,33 DMRG studies found 

that this chiral state can persist also by changing the 
magnetic anisotropy within the XXZ model. ^^’33 in this 
respect, the XXZ model with the only nearest-neighbor 
interactions has not been thoroughly investigated. Some 
recent calculations have pointed out the possibility that 
different magnetic orders are favored for the XY and 
Heisenberg models, based upon order-by-disorder mech¬ 
anisms. This situation should take place for large 
enough spin S', while for S = 1/2 magnetically disor¬ 
dered states should be expected. Indeed, DMRG calcula¬ 
tions have suggested the existence of a spin-liquid phase; 
however, it remains unclear if there is a phase transi¬ 
tion between the XY and the Heisenberg models for 
S = 1/2.3'^’33 In particular, the XY model could have a 
vanishing-small spin gap, which is compatible with a gap¬ 
less quantum spin liquid in the thermodynamic limit. ^3 
Therefore, the spin anisotropy in the nearest-neighbor 
coupling represents a very promising way to unveil the 
nature of the spin-liquid phase of the Heisenberg model. 

In this paper, we consider the following Hamiltonian: 

n = J,y + Sfsp + J.Y ( 1 ) 

{U) iU) 

where {ij) denotes the sum over the nearest-neighbor 
pairs of sites, and = {Sf, Sf, Sf) is spin-1/2 opera¬ 
tor at each site i. In the following, we will set Jxy = 1 
as the energy scale. When Jz = 1 {Jz = 0), Eq. (1) 
reduces to the Heisenberg {XY) model. Here, we focus 
on the region with 0 < < 1, and study the stabil¬ 

ity of different variational wave functions including the 
U{1) and Z 2 spin liquids, as well as the magnetic ordered 
states (the isotropic Heisenberg model with Jz = 1 has 
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FIG. 1: (Color online) (a) The mean-field Ansatzof the U{1) 
Dirac state: black solid (dashed) bonds denote positive (nega¬ 
tive) hoppings (ti). The unit cell is doubled to accommodate 
a magnetic flux $ = tt across hexagons and $ = 0 across tri¬ 
angles. The same hopping amplitudes are also used to define 
the magnetic wave function obtained from Eq. (3). (b) The 
mean-held Ansatz of the Z2[0,7r]/3 state: red solid (dashed) 
lines indicate positive (negative) next-nearest-neighbor hop¬ 
pings (t 2 ) and pairings (A 2 ). 

been thoroughly investigated in previous works^^’^^’^^). 
The main results can be summarized as follow: by ap¬ 
plying an energy optimization, there is no signal for the 
stabilization of a gapped spin liquid; also the inclusion of 
magnetic orders does not improve the U{1) Dirac state. 
Instead, some energy gain can be obtained by including 
a (short-range) Jastrow factor. Then, we construct the 
S = 2 state by exciting 4 spinons in the U{1) Dirac spin 
liquid and calculate the S = 2 spin gap at different val¬ 
ues of Jz. The variational results are compatible with 
the conclusion that the same gapless spin liquid persists 
from Jz = 1 to Jz = 0. 

Method and variational wave funetions. Our varia¬ 
tional wave functions are dehned as 

\^v) = JsVg\^o). ( 2 ) 

where |4^o) is an uncorrelated wave function that is ob¬ 
tained as the ground state of an auxiliary Hamiltonian 
(see below); Vq = Ylii^ ~ is the Gutzwiller pro¬ 

jector, which enforces no double occupation on each site; 
J's = ex.-p{l/2^-j VijSfSj) is the spin Jastrow factor, 
Vij being variational parameters that depend upon the 
distance between sites i and j. We would like to stress 
the fact that such a Jastrow term, which includes the 
^ components of the spin operator, does not break any 
symmetry of the spin Hamiltonian in the easy-plane limit 
{Jz < 1), while it breaks the spin SU (2) symmetry for the 
Heisenberg model {Jz = 1). Here, we consider two cases 
for the auxiliary (non-interacting) Hamiltonian that are 
suitable for magnetic and spin-liquid wave functions. 


Magnetic states are defined from: 

^MAG = ^ Ujcl^^Cj^a- + h ^ Mi • Si, (3) 

(bi),cr i 

where c] ^ (G,cr) creates (destroys) one electron at site 
i with spin a. We find that the best projected state 
within this class of wave functions has non-trivial hop¬ 
ping amplitudes, which define a magnetic fiux 4> = tt 
across hexagons and 4> = 0 across triangles, see Fig. 1 
(they are exactly the ones that define the U{1) Dirac 
spin liquid in Ref. 10). The magnetic order is defined 
by the (variational) parameter h and the vector 
that defines the periodicity; here, we consider copla- 
nar states and restrict in the XY plane, i.e., = 

(cos(r^ • q + sin(r^ • q + r^i), 0) (q is the pitch vector 
and r]i is the phase shift for sites within the same unit 
cell). In the following, we consider two antiferromagnetic 
patterns with q = (0,0) and q = (47r/3,0) (correspond¬ 
ing to the \/3 X \/3 order), see the insets of Fig. 2. With 
Mi in the XY plane, the spin Jastrow factor correctly 
describes the relevant spin fluctuations around the clas¬ 
sical spin state.We would like to emphasize that the 
existence of magnetic long-range order is directly related 
to the presence of a finite parameter h in Eq. (3). 

Instead, spin-liquid wave functions are defined from: 

■^SL = tijcl^Cj,^ + J2[Ajcl^c]^ + h.c.] 

(hi),cr (i,j) 

+ + ( 4 ) 

ia i 

where, in addition to the hopping terms, there is also a 
singlet pairing (A^j = Aji); the on-site paring Aq and the 
chemical potential /i are also considered. It is possible to 
show that many different spin liquids can be constructed, 
depending on the symmetries of Uj and Aij , which may 
have U{1) or Z 2 gauge structure and gapped or gapless 
spinon spectrum.In this paper, we consider two kinds 
of spin liquids, namely the gapless U{1) Dirac state and 



EIG. 2: (Golor online) The hnite size scaling of the varia¬ 
tional parameter /i as a function of the inverse geometrical 
diameter 1 /L at different values of Jz for the magnetic order 
with q = (0, 0) (a) and the a/3 x a/3 order (b). The quadratic 
fitting is used for all cases. The Ansatz for the magnetic order 
with q = (0,0) and a/3 x a/3 are also shown. 
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FIG. 3: (Color online) Energies per site of the U{1) Dirac 
spin liquid (green diamonds) compared with the one of the 
magnetic state with q = (0,0) (red circles) for different values 
of Jz- The U{1) Dirac state with spin Jastrow factor is also 
reported (blue squares), which corresponds to the magnetic 
state with a vanishing variational parameter h. 


the gapped Z2[0,7r]/3 state, as shown in Fig. 1. We em¬ 
phasize that, for < 1, a spin Jastrow factor can be 
also included, since the spin SU{2) symmetry is explic¬ 
itly broken by the XXZ Hamiltonian (1). 

In general, suitable boundary conditions in the auxil¬ 
iary Hamiltonians are chosen, in order to have a unique 
mean-field ground state order to fulfill the con¬ 

straint of one electron per site (imposed by the Gutzwiller 
projector) and to take into account the spin Jastrow fac¬ 
tor, a Monte Carlo sampling is needed. Along the Markov 
chain that define the numerical simulation, only configu¬ 
rations belonging to the physical Hilbert space are pro¬ 
posed, so that the Gutzwiller projector is exactly im¬ 
plemented. To optimize the variational parameters in 
both the auxiliary Hamiltonians and the spin Jastrow 
factor, we use the stochastic reconfiguration (SR) opti¬ 
mization method to find the energetically favored state 
in variational Monte Carlo scheme.The SR optimiza¬ 
tion method allows us to perform the optimization with 
many variational parameters, and to obtain an extremely 
accurate determination of variational parameters. 

Results. We performed the variational Monte Carlo 
calculations on toric clusters with Lx Lx3 sites and pe¬ 
riodic boundary conditions. Let us start with magnetic 
states. In Fig. 2, we show the size scaling of the mag¬ 
netic order parameter h of Eq. (3) for both the states 
with q = (0, 0) and \/3 x ^/3 order. The latter one is 
not frustrated by boundary conditions only when L is a 
multiple of 3. First of all, we find that a finite magnetic 
parameter h of Eq.(3) can be stabilized on finite clus¬ 
ters, and the two states have essentially the same energy 
within 10~^Jxy for all the cases analyzed here. How¬ 
ever, the most important outcome is that ^ 0 in the 
thermodynamic limit (for both antiferromagnetic states), 
indicating that no magnetic order can be stabilized in the 
XXZ model. Therefore, the auxiliary Hamiltonian from 


which the variational state is constructed reduces to the 
U{1) Dirac state of Ref. 10. Nevertheless, the param¬ 
eters Vij of the spin Jastrow factor remain finite, with 
sizable values at short-range distances. This fact gives 
a non-negligible energy gain with respect to the U{1) 
Dirac state, especially close to the XY limit, see Fig. 3. 
We stress that the presence of the spin Jastrow factor 
is just relevant to improve short-range observables, such 
as the energy. We mention that the spin Jastrow factor 
gives a small energy gain of about 10“^ also in the 
Heisenberg case, see Fig. 3 (however, in this case, the 
spin Jastrow factor spoils the spin SU{2) invariance of 
the Dirac spin liquid). 

We now move to study the possible stabilization of 
Z 2 spin liquid states. According to the classification of 
Ref. 42, there is only one Z 2 gapped spin liquid that is 
directly connected to the U{1) Dirac spin liquid, the so- 
called ^2[0, tt]^ state. In Ref. 11, the authors have shown 
that this gapped spin liquid cannot be stabilized in the 
Heisenberg model with = 1. Here, we would like to 
extend the analysis to the case of the XXZ model. The 
variational states is constructed from Eq. (4), and the 
non-interacting wave function that defines the Z2[0,7r]/3 
state includes the nearest-neighbor ti, the next-nearest- 
neighbor t 2 hoppings, the next-nearest-neighbor paring 
A 2 , which is responsible for the breaking from U{1) to Z 2 
symmetry, a chemical potential /i, and the on-site pair¬ 
ing Aq, see Eig. 1. In the following, we do not consider 
the spin Jastrow factor, which may improve the energy 
but does not change the optimization of the variational 
parameter A 2 . The optimization is shown in Eig. 4 for 
Jz = 0, 0.2, and 0.6 on the L = 16 cluster. The result 
is that, as for the Heisenberg case, both A 2 and Aq (not 
shown here) go to zero for all the values of Jz consid¬ 
ered, even on finite clusters. The vanishing Z 2 parame¬ 
ters indicates that, similar to what has been found in the 



FIG. 4: (Golor online) The variational parameter A 2 for 
the Z2[0,7r]/3 spin liquid for Jz = 0, 0.2, and 0.6 on L = 16 
cluster. The variational energy per site as function of the SR 
iterations is shown in the inset. 
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FIG. 5: (Color online) The finite size scaling of the S = 2 
spin gap as a function of 1/L for different values of Jz- The 
results for Jz — 1 are from Ref. 26. 

Heisenberg model,the gapped Z2[0,7r]/3 spin liquid is 
not stable in the XXZ model. 

In summary, we obtain that the best variational state 
of the form (2) that can be constructed from Eqs. (3) 
and (4) is the U{1) Dirac spin liquid, supplemented by a 
short-range spin Jastrow factor.^^ 

In the following, we compute the S = 2 spin gap for 
the U{1) Dirac state. The S = 2 state is constructed by 
changing boundary conditions, in order to have 4 spinous 
in an 8-fold degenerate single-particle level at the chem¬ 
ical potential; a unique mean-field state is then obtained 
by taking all these spinons with the same spin. This 
S = 2 state can be written in terms of a single deter¬ 
minant, which is particularly easy to be treated within 
our Monte Carlo sampling. In the following, we do not 
consider the spin Jastrow factor, since its inclusion does 
not modify the qualitative picture.Similarly to what 
has been done on the Heisenberg model,we obtain the 
S = 2 spin gap by computing separately the energies of 
the S = 0 and S = 2 states. In Fig. 5, we report the 
results for Jz = 0, 0.2, and 0.6 (the case with Jz = 1 
from Ref. 26 are also reported for comparison). First of 
all, we remark that, for each cluster size, the spin gap de¬ 
creases by decreasing the value of Jz , indicating that the 
anisotropy in the spin super-exchange tends to close the 
finite-size gap. This result is in agreement with DMRG 


calculations in Ref. 38. Most importantly, the finite-size 
scaling with L up to 20 clearly indicates a vanishing spin 
gap for all values of Jz considered here. Therefore, our 
analysis based upon Gutzwiller-projected states suggests 
that the same U{1) Dirac state with gapless spinon exci¬ 
tations can be stabilized from = 1 to = 0. 

Conclusions. In summary, we investigated the XXZ 
model on the kagome lattice by using the variational 
Monte Carlo technique with the Gutzwiller projected 
fermionic states. We have studied different variational 
wave functions describing either magnetic states or spin- 
liquid phases. As previously obtained in the Heisenberg 
model,the gapped Z2[0,7r]/3 spin liquid cannot be sta¬ 
bilized for Jz < indicating a remarkable stability of 
the gapless U{1) Dirac state. Moreover, the considera¬ 
tion of magnetic order does not give any energy gain in 
the thermodynamic limit, but a considerable energy gain 
can be obtained by the spin Jastrow factor. The S = 2 
spin gap, on any finite-size clusters, decreases with de¬ 
creasing the value of L^, indicating that the best place 
to find a gapless spin liquid is most probably close to the 
XY limit. This outcome agrees with a recent DMRG 
study^^, which suggested that a critical state can be sta¬ 
bilized near the XY kagome model. In addition, we do 
not find any evidence for possible dimer states, as also 
suggested by DMRG calculations. 

Finally, we would like to mention that the application 
of few Lanczos steps to the 17(1) Dirac spin liquid, as al¬ 
ready done in recent works for the Heisenberg model, 
does not alter the results on the S = 2 gap, although the 
large statistical errors for L = 8 do not allow us to obtain 
as neat conclusions as in the Heisenberg model.In fact, 
even though the U{1) Dirac spin liquid remains the best 
variational state within the class of fermionic states that 
have been analyzed here, its accuracy slightly deterio¬ 
rates when decreasing the value of Jz^ which makes the 
zero-variance extrapolation harder than for the Heisen¬ 
berg case. 
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FIG. 6: (Color online) The parameters Vij as function of the 
distance \i — j\ for the L — 12 lattice at Jz = 0, 0.2, 0.6, and 
1.0. The black lines are guides for eye. 

Supplemental Material 

The spin Jastrow factor- The inclusion of the spin Jas- 
trow factor Js = ex.-p{l/2^-j VijSf Sj) gives rise to a 
considerable improvement of the U{1) Dirac state in the 
whole region 0 < < 1, as we discussed in the main 

text. In Fig. 6, we show the optimized parameters Vij 
as function of the distance \Ri — Rj\ between the sites 
i and j. All parameters are positive and decay expo¬ 
nentially with distance, indicating that the spin Jastrow 
factor is short range. Moreover, the rate of decay seems 
to increase for \Ri — Rj\ > 9. 

The Lanczos steps- In the main part of the paper, we 
concentrate on variational wave functions as defined by 
Eq. (2); moreover, in few cases, we also improve them by 
applying a number p of Lanczos steps: 

i^p)= + (5) 

where am are p additional variational parameters. 
Clearly, whenever is not orthogonal to the exact 

ground state, |^p) converges to it for large p. Unfortu¬ 
nately, on large sizes, only few steps can be efficiently 
afforded: here, we consider the case with p = 1 and 
p = 2 {p = 9 corresponds to the original variational 
wave function). Furthermore, an estimate of the exact 
energy may be obtained by the variance extrapolation. 
Indeed, for a systematically convergent sequence of states 


with energy Ep and variance cr^, it is easy to prove 
that Ep ^ Eex + const x where Ep = {^p\H\il/p)/N 
and are the energy 

and variance per site, respectively. Therefore, the ex¬ 
act energy Eq^ may be extracted by fitting Ep vs cr^, for 
p = 0,1, and 2. 

Few Lanczos steps may be applied to the variational 
states with S = 9 and S' = 2, as described in Eq. (5), 
allowing a zero-variance extrapolation of the energies. 
Since this procedure is quite computational demanding, 
we only consider the XY limit. Let us start by dis¬ 
cussing the results on a small system with L = 4. Eor 
S = 0, the p = 2 state has E = —0.299204(2), while, 
by performing the zero-variance extrapolation, the esti¬ 
mated energy is E = —0.30045(1), which is quite close to 
the DMRG result E = —0.301228 on long cylinder (see 
Eig. 7a). Eor the S = 2 excitation, the p = 2 energy is 
E = —0.295072(2) and the zero-variance extrapolation 
gives E = —0.29668(3), again very close to the DMRG 
result of E = —0.29'744 (see Eig. 7b). Remarkably, the 
extrapolated gap that we obtain agrees with the DMRG 
one, indicating that there is an almost exact cancellation 
error between the S = 9 and S = 2 energies. 

All energies with p = 0, 1, and 2 Lanczos steps and 
their zero-variance extrapolation on L = 4, 6, and 8 clus¬ 
ters are reported in Table I and Eig. 7(a and b). Com¬ 
pared to the Heisenberg model {Jz = 1),^^ the S = 2 
spin gaps on finite size clusters (L = 4, 6, and 8) are 
smaller (Eig. 7(c)). This is consistent with the DMRG 
calculations.^^ However, the variance in the XY model is 
almost twice time larger than the one in the Heisenberg 
model for each Lanczos step. This fact indicates that the 
1/(1) Dirac state is less accurate to describe the ground 
state when decreasing Jz. Nevertheless, the thermody¬ 
namic extrapolation of the S = 2 gap is still possible, see 
Eig. 7(c). Here, the large error bar of the L = S cluster is 
entirely due to the large variances of the S = 9 and S = 2 
states, which makes a rather imprecise extrapolation of 
zero variance. By performing a fit of the three sizes with 
L = 4, 6 and 8, which takes into account their error bar, 
we obtain a vanishing S = 2 spin gap in the thermody¬ 
namic limit. Taking into account all the statistical errors 
of the fitting procedure, the largest possible value for the 
thermodynamic gap is about 0.05. 

In the main text we have shown that, the spin Jastrow 
factor improves the ground-state energy of the pure 1/(1) 
Dirac state. We have also performed the Lanczos steps 
on this wave function for the XY model, and obtained 
the smaller variance with p = 9 and 1. However, the 
p = 2 calculation is unstable: with a small change in the 
p = 2 Lanczos parameters, the variance may have large 
variations, while the energy does not change much. This 
fact may indicate that there are many low-lying states 
with competing energy for the XY model. Nevertheless, 
on the L = 4 cluster, by performing a linear extrapolation 
with p = 9 and 1 results, we get the DMRG energies 
within one error bar for both the ground state and S = 2 
excitation. 
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TABLE I: Energies of the U{1) Dirac spin liquid (columns 2-5) and its S' = 2 excitation (columns 6-9), with p = 0, 1, and 2 
Lanczos steps on different clusters for the spin-1/2 XY model. The estimated energies of the iS = 0 and S = 2 states by using 
the zero-variance extrapolation marked in bold. 



p = 0 p=l p = 2 S = 0 

p = 0 p=l p = 2 S = 2 

II II II 

00 Oi 

-0.2862336(7) -0.2966459(4) -0.299204(2) -0.30045(1) 
-0.2858440(5) -0.2947757(7) -0.297948(1) -0.30024(5) 
-0.2857821(6) -0.2934697(6) -0.296803(3) -0.3002(1) 

-0.2816427(7) -0.2923336(5) -0.295072(2) -0.29668(3) 
-0.2843146(5) -0.2934122(6) -0.296711(2) -0.29921(5) 
-0.2850986(5) -0.292854(1) -0.296242(3) -0.2998(2) 




EIG. 7: (Color online) Energies per site for the S = 0 ground state (a) and S = 2 excitation (b) versus the variance for Jz = 0. 
The results with p = 0, 1, and 2 are reported for L = 4, 6, and 8. The variance extrapolated results are shown. The DMRG 
results on long cylinder with L = 4 are also reported, (c) The S = 2 spin gap with p = 0 and extrapolation as a function of 
the inverse geometrical diameter (l/L) at Jz = 0. The U{1) Dirac state with and without Jastrow factor are both considered. 
The results for Jz — 1 are from Ref. 26. On L = 4 cluster, the results for the spin gap obtained by different wave functions are 
the same and equal to the one obtained by DMRG. Eor the XY model, the upper bound of the S = 2 spin gap is 0.05 in the 
thermodynamic limit: this is entirely due to the large statistical error on L = 8. 


















